"""
Peak detection
"""
import operator
import atlas, numpy
from atlas import util, bitvector
from itertools import izip, count, islice, imap, ifilter, tee, repeat

@util.gc_off
def make_list( size, default=0, pad=0, ntype=numpy.int8):
    
    #
    #uncomment the line below to drastically reduce memory use
    #at the cost of increasing runtime with an order of magnitude
    #

    #return bitvector.BitVector( size=size+pad, default=default )

    if default==0:
        return numpy.zeros( size+pad, ntype )
    else:
        return numpy.ones( size+pad, ntype )
        

def detect_peaks( data ):
    """
    Detects peaks (local maxima) from an iterator that generates 
    float values. Input data is a list of tuples containing 
    the index and the value at the index.
    
    Returns a generator yielding tuples that corresponds 
    to the peak index and peak value.
    
    >>> values  = [ 0.0, 1.0, 2.5, 1.0, 3.5, 1.0, 0.0, 0.0, 10.5, 2.0, 1.0, 0.0 ]
    >>> data  = izip( count(), values )
    >>> peaks = detect_peaks( data )
    >>> peaks = list(peaks)
    >>> peaks
    [(2, 2.5), (4, 3.5), (8, 10.5)]
    >>> select_peaks( peaks, width=1)
    [(2, 2.5), (4, 3.5), (8, 10.5)]
    >>> select_peaks( peaks, width=2)
    [(4, 3.5), (8, 10.5)]
    """

    # a more functional approach using less memory
    def is_local_maxima( row ):
        "Selects local maxima"
        return row[0][1] < row[1][1] >= row[2][1]

    teedata   = tee(data, 3)
    mapped    = imap( islice, teedata , count(), repeat(None) )
    zipped    = izip( *mapped )
    filtered  = ifilter( is_local_maxima, zipped )
    for left, mid, right in filtered:
        yield mid

def select_peaks( peaks, width, level=0):
    """
    Selects successive non-overlapping peaks 
    with a given exclusion zone.

    Takes as input a list of (index, value) tuples.

    Returns a list of tuples = ( midpoint, maxima )
    """
     
    # order by peak height
    work  = [ (y, x) for x, y in peaks if y >= level and x > width ]
    work.sort()
    work.reverse()

    # sanity check
    if not work:
        return []

    #largest index        
    xmax = max(work, key=operator.itemgetter(1) )[1]
    
    # selected peaks
    selected = []

    # keeps track of empty regions
    #empty = bitvector.BitVector( size=xmax+width+1, default=1 )
    
    empty = make_list( size=xmax, pad=width+1, default=1 )
    
    half = width/2
    for peaky, peakx in work:
        if empty[peakx]:
            left  = peakx - width
            right = peakx + width
           
            # store into output data
            selected.append( ( peakx, peaky ) )
            
            # block the region
            empty[left:right] = numpy.zeros (right - left)

    selected.sort()
    return selected

def _test():
    """
    Main testrunnner
    """
    import doctest
    doctest.testmod( optionflags=doctest.ELLIPSIS + doctest.NORMALIZE_WHITESPACE )
    
if __name__ == "__main__":
    _test()

    values  = [ 0.0, 1.0, 2.5, 1.0, 3.5, 1.0, 0.0, 0.0, 10.5, 2.0, 1.0, 0.0 ]
    data  = izip( count(), values )
    peaks = detect_peaks( data )
    peaks = list(peaks)
    
    print select_peaks( peaks, width=2)